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Abstract 

We  develop  new  quadrature  rules  for  Isogeometric  Analysis  based  on  the  solution  of  a  local  nonlin¬ 
ear  problem.  A  simple  and  robust  algorithm  is  developed  to  determine  the  rules  which  are  exact 
for  important  B-Spline  spaces  of  uniform  and  geometrically  stretched  knot  spacings.  We  consider 
both  periodic  and  open  knot  vector  configurations  and  illustrate  the  efficiency  of  the  rules  on  se¬ 
lected  boundary  value  problems.  We  find  that  the  rules  are  almost  optimally  efficient,  but  much 
easier  to  obtain  than  optimal  rules,  which  require  the  solution  of  global  nonlinear  problems  that 
are  often  ill-posed. 

Keywords:  Numerical  integration,  Isogeometric  analysis,  NURBS,  B-splines 


1.  Introduction 

Isogeometric  analysis  (IGA)  is  a  recently  proposed  computational  technique  [19]  for  the  solution 
of  boundary  value  problems,  based  on  the  idea  of  using  the  same  functions  adopted  in  Computer 
Aided  Design  (CAD)  not  only  to  describe  the  domain  geometry,  but  also  to  build  the  numerical 
approximation  of  the  problem  solution. 

In  CAD  a  common  choice  for  the  representation  of  complex  geometries  is  Non-Uniform  Rational 
B-Splines  (NURBS).  Thus,  a  particular  case  of  isogeometric  methods  is  to  represent  the  geometry, 
project  the  data,  and  find  the  solution  in  the  space  of  NURBS.  For  the  interested  reader,  we  recall 
that  IGA  has  been  summarized  in  a  recent  book  [13],  studied  in  a  number  of  contributions  (e.g., 
[2,  4,  7,  14,  16,  17,  20]),  and  that  it  is  having  a  progressively  growing  impact  on  fields  as  diverse 
as  fluid  dynamics  [5,  6,  8,  18],  structural  mechanics  [1,  3,  15,  22,  24],  and  electromagnetics  [10,  9]. 

The  reason  for  such  extensive  attention  is  mainly  due  to  some  important  advantages  of  IGA 
methodology  over  more  classical  Finite  Element  Analysis  (FEA).  An  example  is  the  fact  that  IGA 
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numerical  solutions  can  possess  high  regularity  across  mesh  elements.  This  is  an  important  feature, 
leading  to  a  higher  accuracy  per  degree-of- freedom  (see  [7,  17]),  to  improved  spectrum  properties  of 
the  discrete  operators  (see  [20]),  and  the  possibility  of  constructing  discretizations  able  to  preserve 
fundamental  structures  of  the  continuum  differential  operators  (such  as  De  Rham  diagrams,  see 

[9])- 

However,  within  IGA  approaches  a  remaining  issue  is  the  design  of  efficient  quadrature  rules 
taking  advantage  of  high  inter-element  regularity.  In  fact,  the  technique  adopted  heretofore  in  the 
literature  is  to  use  Gauss  quadrature  on  each  element,  a  choice  far  from  being  optimal. 

The  aim  of  the  present  paper  is  to  design  an  efficient  quadrature  strategy,  in  the  sense  that 
it  should  result  in  a  considerable  savings  in  terms  of  computational  effort  compared  to  classical 
Gauss  rules,  while  maintaining  the  classical  element-by-element  assembly  procedure. 

The  construction  of  quadrature  rules  for  IGA  was  initially  considered  in  [21].  The  rule  proposed 
therein  is  optimal  because  it  exactly  integrates  B-spline  basis  functions  with  the  minimum  number 
of  function  evaluations;  however,  such  an  optimal  quadrature  rule  is  obtained  as  a  solution  of  a 
global,  non-linear  and,  in  general,  ill-conditioned  system  of  equations,  which  is  therefore  difficult 
to  solve  for  high  polynomial  degrees  and  numbers  of  elements.  Moreover,  such  a  quadrature  rule 
has  to  be  recomputed  when  the  number  of  elements  changes.  As  a  consequence,  in  [21]  the  authors 
suggest  utilizing  the  rule  in  a  non-optimal  way,  considering  the  function  regularity  only  at  the  level 
of  small  groups  of  elements  (so-called  macro-elements). 

The  present  paper  takes  advantage  of  the  translation-invariance  of  the  basis  functions  defined 
in  the  domain  interior,  leading  to  a  quadrature  rule  obtained  solving  a  local,  non-linear  system 
of  equations,  which  does  not  have  to  be  recomputed  when  the  number  of  elements  varies.  As 
a  consequence,  the  new  rule  competes  well  with  the  optimal  one  proposed  in  [21]  in  terms  of 
efficiency,  but  is  much  simpler  to  construct.  However,  special  considerations  needs  to  be  given  to 
elements  on  the  boundary. 

The  rule  is  initially  designed  for  the  case  of  a  uniform  discretization  of  elements  in  the  para¬ 
metric  domain,  but  is  then  extended  to  the  case  of  geometrically  refined  meshes. 

The  paper  is  organized  as  follows.  Section  2  describes  some  preliminaries  on  IGA  and  quadra¬ 
ture.  Section  3  discusses  one- dimensional  B-splines,  and  presents  the  new  quadrature  rule,  first 
for  a  periodic  uniform  knot  vector  and  then  for  an  open  uniform  knot  vector,  detailing  the  algo¬ 
rithm  and  discussing  optimality.  Section  4  exploits  the  new  quadrature  rule  to  numerically  solve 


2 


boundary  value  problems,  including  one  whose  solution  possesses  sharp  boundary  layers.  Section 
5  gives  some  final  comments. 


2.  Preliminaries  on  IGA  and  quadrature 


For  the  sake  of  simplicity,  the  following  presentation  focuses  on  a  two-dimensional  setting,  which 
can  be  easily  converted  to  situations  with  a  different  dimensionality.  Accordingly,  introducing  a 
parametric  domain  12  =  [0,  k]  x  [0,  A;],  the  solution  of  classical  engineering  problems  requires  the 
calculation  of  integrals  such  as 


Jn 

(2.1) 

(2.2) 

[  viM*)i2j(0*(0<& 

Jh 

(2.3) 

where:  i  and  j  are  two-dimensional  multi-indices;  and  /?j  are  NURBS  basis  functions;  <f>  is 
a  factor  taking  into  account  the  coefficients  of  the  investigated  partial  differential  equation  and 
the  problem  geometry,  i.e.,  the  Jacobian  of  the  geometry  map  and,  possibly,  the  derivatives  of  its 
inverse  [13].  Typically,  (2.1)  emanates  from  a  mass  term,  (2.2)  from  a  stiffness  term,  while  (2.3) 
corresponds  to  a  linear  advection  term,  classically  encountered  in  fluid  dynamics.  Of  course,  one 
can  consider  also  other  terms,  which  could  be  included  in  the  forthcoming  discussion  without  any 
difficulty. 

As  commented  in  [21,  Section  3.2],  within  isogeometric  analysis  it  is  a  common  practice  to 
compute  terms  like  (2.1)-(2.3)  adopting  quadrature  rules  able  to  exactly  integrate  the  spline  basis 
functions  that  generate  the  NURBS.  Furthermore,  in  the  case  of  tensor-product  discretizations,  the 
multi-dimensional  quadrature  rule  is  obtained  from  the  tensor-product  of  suitable  one-dimensional 
rules. 

The  one- dimensional  integrals  of  interest  are  of  the  following  kind 

k 

(pi(x)(j)j(x)  dx, 


<l>i(x)<f>j(x)dx, 


(j)'i(x)(j)j(x)dx, 


(2.4) 
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where  4>i  and  (f)j  belong  to  a  proper  space  of  univariate  spline  functions.  Assuming  the  use  of  the 
approximation  space  S ft,  i.e.,  the  space  of  spline  functions  of  degree  p  and  regularity  r,  examining 
Equation  (2.4),  we  may  conclude  that  for  the  design  of  an  exact  quadrature  rule  the  space  of 
integrand  functions  to  be  considered  is  S^_v 

Considering  as  an  example  the  interesting  case  of  shape  functions  with  maximum  regularity, 
i.e.,  G  Sp_ll  we  have: 


e  Sp_i, 

m  6  S?_f,  (2.5) 

•VAj  6  s?_-2\ 

and,  therefore,  »S'|P2  is  the  space  that  contains  all  the  integrands. 

To  generalize  the  discussion,  in  the  following  we  indicate  the  space  of  the  functions  to  be 
integrated  as  S'™,  where  m  denotes  the  degree  of  the  functions  to  be  integrated  and  q  denotes  the 
corresponding  regularity;  clearly,  the  quantities  m  and  q  come  from  the  original  problem  of  interest 
(i.e.,  from  Equations  (2.4)  and  (2.5)),  hence  they  should  satisfy  relation 

m  =  2p  ,  q  <  p  —  1, 


from  which  it  can  also  be  derived 

M  -1- 

where  [s]  is  the  smallest  integer  l  such  that  s  <  l.  Even  though  in  this  context  m 
m  £  N  for  the  sake  of  generality. 


(2.6) 


2 p,  we  allow 


3.  Quadrature  of  univariate  splines 


Consistent  with  the  previous  discussion,  we  now  aim  to  construct  a  quadrature  rule  able  to 
exactly  evaluate  an  integral  in  the  form 

ipi(x)dx 

for  the  case  of  a  spline  function  ?/y  G  S^,  where  m  and  q  are  related  by  (2.6). 

To  this  end  we  recall  that,  given  a  generic  function  /,  an  n-point  quadrature  rule  is  a  choice  of 
n  ordered  points  and  weights  ( Xi,Wi )  such  that 

b  n 

f(x)  dx  ~  J2wJ(Xi). 

1=1 
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A  quadrature  rule  is  said  to  be  exact  on  the  family  of  functions  {fj}j=i,...,M  if 

b  n 

fj(x )  dx  =  Y^  Wifj(xi)  Vj  =  1, ,  M.  (3.1) 

i—  1 

Condition  (3.1)  represents  M  non-linear  equations  in  the  2 n  unknowns  ( Xi,Wi ).  Classically,  exact¬ 
ness  is  required  on  polynomials  np  to  a  fixed  degree;  in  particular,  the  n-point  Gauss  rule  is  the 
(only)  n-point  rule  that  is  exact  for  polynomials  np  to  degree  2 n  —  1. 

Recently,  exactness  on  general  families  of  functions  has  been  also  explored.  For  example,  in 
[23,  12]  a  quadrature  rule  that  integrates  exactly  M  independent  functions  {/],}  using 
is  introduced  and  named  Generalized  Gaussian  with  respect  to  the  functions  {fj}. 

With  the  previously  addressed  general  consideration,  to  properly  approach  the  problem  of  our 
interest,  in  the  following  we  provide  some  details  on  one- dimensional  B-splines  (Section  3.1).  We 
then  approach  the  quadrature  problem,  for  a  periodic  uniform  knot  vector  (Section  3.2)  and  for 
an  open  uniform  knot  vector  (Section  3.3),  presenting  also  an  algorithmic  description  (Section  3.4) 
and  a  discussion  on  the  optimality  (Section  3.5)  for  the  newly  proposed  quadrature  rule. 

3.1.  B-splines  in  one  space  dimension 

As  a  one- dimensional  parametric  domain,  referred  to  as  a  patch ,  we  assume  the  interval  [0,  k] 
with  a  uniform  subdivision  into  unitary  elements  \j  —  1,  j],  where  j  —  1, ...  ,k. 

A  spline  function  /  G  S™  is  a  polynomial  of  degree  m  in  each  element  [j  —  1,  j]  with  /  € 
C9([0,A;]).  Accordingly,  to  determine  a  function  /  e  S™,  we  need  to  assign  m  +  1  polynomial 
coefficients  on  each  of  the  k  elements  with  q  +  1  continuity  requirements  on  the  k  —  1  internal 
points;  therefore,  indicating  by  N  the  dimension  of  the  space  5™,  the  following  relation  holds: 

N  =  k{m  +  1)  —  (k  —  l)(q  +  1)  =  k{m  —  q)  +  q  +  1. 

In  general,  q  <  m  —  1  is  allowed,  even  though  in  this  paper  we  further  restrict  to  condition  (2.6). 
We  also  assume  q  >  —  1,  where  q  —  —  1  refers  to  the  discontinuous  case,  while  q  =  0  refers  to  the 
case  of  functions  continuous  on  the  whole  domain  [0,  k\  and  piecewise  polynomial  on  each  single 
element,  i.e.,  the  classical  FEA  approximation. 

Following  a  standard  procedure  in  the  B-spline  literature  [25],  to  construct  the  basis  functions 
we  define  a  knot  vector  in  the  form 


M 

Y 


points 


{si;  ^2;  •  •  •  ;  £/V+m+l} 
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In  particular,  as  a  first  case  we  consider  an  open  uniform  knot  vector  H,  defined  as 


s  =  {(WA  w,  •  •  .,fc-l,...,fc-lfc,...,fc}  (3.2) 

ra+1  times  m—q  times  m—q  times  m+1  times 

where  the  hrst  and  last  knots,  i.e.,  0  and  k,  have  multiplicity  (m  +  1),  and  the  other  knots  have 
multiplicity  (m  —  q).  Starting  from  5,  we  can  construct  the  basis  functions  for  the 

space  5+  using  the  following  recursive  formula 


(+  =  7 - 1  0)  +  7 — $j+i (x) >  1  ~  1 . 


with  the  initial  condition 


+°+)  = 


£j+^+i  —  £j+i 

1  if  ij  <  x  <  0+1 
0  otherwise 


(3.3) 


and  setting 


OjOO  = 


We  recall  the  convention  that,  when  a  denominator  is  zero  in  (3.3),  the  corresponding  quotient  is 
set  to  zero. 

Similarly,  we  can  introduce  a  periodic  uniform  knot  vector  H,  defined  as 


s  =  {01_^0,  , . . .  ,fc  — I,---,*:  — 1,  (3.4) 

m—q  times  m—q  times  m—q  times  m—q  times 

and,  adopting  the  recursive  scheme  indicated  in  (3.3),  we  can  construct  the  basis  functions  asso¬ 
ciated  to  S  and  spanning  a  space  denoted  as  S™.  We  observe  that  S™  C  S™  and,  in  particular, 
that  =  span  {i^j}j=q+2i  N_q_v 

We  recall  that  with  this  construction  each  basis  function  is  non-negative  and 


supp ipj  =  [&,fj+m+i] ; 

therefore,  since  for  the  specific  cases  of  interest  we  have  an  m  —  q  multiplicity  of  internal  knots 
(e.g.,  Equations  (3.2)  and  (3.4))  as  well  as  restriction  (2.6)  on  m  and  q,  each  B-spline  basis  function 
ij)j  can  have  support  in  at  most  two  elements. 

Accordingly,  we  introduce  a  new  notation,  to  distinguish  between  basis  functions  having  support 
only  in  one  element  and  basis  functions  having  support  on  two  elements.  In  particular,  we  indicate 
with  the  superscript  (i)  the  basis  functions  that  have  support  included  in  the  element  [i  —  1,  i]  and 
with  the  superscript  (i,i  +  1)  those  with  support  in  [i  —  1,  i  +  1], 
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Moreover,  we  observe  that  the  basis  functions  ifj  G  S™  with  support  only  on  one  element  can 
be  further  distinguished  into  a  set  of  boundary  functions,  in  the  following  also  indicated  with  a 
superposed  bar,  which  are  present  only  in  the  first  and  the  last  element,  as  a  consequence  of  the 
fact  that  the  first  and  the  last  knot  have  multiplicity  m  +  1,  and  a  set  of  bubble  functions,  which 
are  instead  present  in  all  the  elements,  i.e.,  both  in  the  boundary  as  well  as  in  the  internal  ones. 
Finally,  all  the  functions  with  support  on  two  elements  are  referred  to  as  transmission  functions. 

The  whole  set  of  basis  functions  for  the  space  S'™  can  be  rearranged  distinguishing  between 
three  groups  of  functions: 

•  2 (g  +  1)  boundary  functions  ,  1  <  j  <  q  +  1,  i  —  1,  k; 

•  k(m  —  2q  —  1)  bubble  functions  ,  1  <  j  <  m  —  2q  —  1 ,  1  <  i  <  k; 

•  (k  —  l)(g  +  1)  transmission  functions  VqM+ '  —  1- 

As  an  example,  for  the  case  m  =  10,  k  —  3  with  q  =  3,  Figure  1  reports  all  the  basis  functions 
generated  from  the  open  uniform  knot  vector,  while  Figure  2  reports  the  same  basis  functions, 
distinguishing  between  boundary,  bubble  and  transmission  functions.  We  observe  that 

•  the  q  +  1  boundary  functions  on  the  first  element  are  related  to  the  q  +  1  boundary  functions 

on  the  last  element  by  the  following:  =  ^^2-j(k  —  x )  (we  wi^  saY  that  they  are 

specular) ; 

•  the  bubble  functions  are  such  that  V*i,*2  G  {1, the  following  translation  property 
holds:  ifjll\x)  =  i>j  2\x  -  (ii  -  i2)); 

•  the  transmission  functions  are  such  that  VAA2  £  {l,---,k  —  1}  the  following  translation 
property  holds:  if^lM+l\x)  =  ipjl2,l2+1\x  —  (i\  —  i2)). 

Recalling  that  S™  C  S™ ,  we  observe  that  the  space  Sff  is  spanned  by: 

•  k(m  —  2q  —  1)  bubble  functions  ,  1  <  j  <  m  —  2q  —  1 ,  1  <  i  <  k\ 

•  (k  —  l)(g  +  1)  transmission  functions  ,  l<j<q  +  l,l<i<k  —  1. 

At  this  stage  it  is  fundamental  to  observe  that  all  functions  in  S™  are  translation  invariant,  and 
this  property  will  be  used  in  the  construction  of  an  optimal  rule  for  exact  quadrature  in  Section 
3.2.  Since  S™  and  C  S'™  differ  only  by  the  presence  of  the  boundary  functions  in  S'™,  in  Section 
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Figure  1:  Basis  functions  generated  by  the  open  uniform  knot  vectors  for  the  case  m  =  10,  k  =  4,  q  =  3. 


3.3  we  will  reach  our  final  goal  of  developing  an  exact  quadrature  for  S™,  enriching  the  rule  for 
S™  in  order  to  integrate  the  additional  boundary  functions. 


3.2.  Quadrature  on  a  periodic  uniform  knot  vector 

The  first  problem  we  want  to  address  is  the  construction  of  an  exact  quadrature  rule  for 
fjj  £  S™.  As  we  have  seen  in  the  previous  section,  these  functions  can  be  grouped  in  bubble 
and  transition  functions  that  are  repeated  up  to  translation  on  each  element  or  couple  of  elements, 
respectively.  Because  of  the  translation  properties  of  the  functions  to  be  integrated,  we  also  assume 
that  the  quadrature  rule  is  the  same  in  all  elements.  That  is,  for  i  £  {1, . . . ,  k}  and  1  =  1,...,  n1 , 
if  x'l’1  and  wj’1  are  the  quadrature  points  in  [i  —  1  ,i\  and  quadrature  weights,  respectively,  we  take 

XY  =  xi  +  (*  “  !),  and  wi’1  =  wi , 


with  xj  £  [0,1].  The  number  of  quadrature  points  per  element  is  denoted  by  n1 ,  and  will  be 
determined  later.  Then,  the  problem  of  the  exact  integration  of  the  bubble  functions  fjj 

nI  M, 

( x\ ,l)  =  /  i/jj'\x)  dx,  1  <  j  <  m  —  2q  —  1, 

i=i  ^ *-1 

and  exact  integration  of  the  transition  functions  (for  i  £  {1 ,...  ,k  —  1}) 

n1  n1 

i=i 


i=i 


I  i— 1 


b.*+1)(a;)  fjx^  1  <  j  <  q  +  1 


can  be  summarized  as  follows 


1 


Figure  2:  Basis  functions  generated  by  the  uniform  open  knot  vector  for  the  case  m  =  10,  k  =  4, 
between  boundary  functions,  bubble  functions  and  transmission  functions,  respectively. 


q  =  3,  distinguishing 
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Problem  3.1  (Quadrature  for  Periodic  Problem).  For  a  given  i  G  {1, ...  k  —  1},  find  quadrature 
points  xj  G  [0,1],  and  weights  w{  >  0  ,  l  —  1, . . . ,  n1  such  that 


wiy\x\  +  (i  —  1))  =  /  y\x)  dx ,  1  <  j  <  m  —  2q  1 

i=i  i~1 

n1 

E 

i=i 


w, 


i+ 1 


(xj  +  (i-  1))  +-0  i  ’(xj  +  i) 


(3.5) 


1  <  j  <  q  +  I- 


Vu-v+1)  (x)  dx 


'  i— 1 


Since  (3.5)  is  a  system  of  m  —  q  equations,  we  look  for  n1  = 
weights. 


(m  —  g) 


□ 


quadrature  points  and 


In  particular,  for  the  case  (m  —  g)  even,  we  have  a  square  non-linear  system  of  equations  and, 
as  shown  in  the  numerical  simulations  in  Section  4,  this  leads  to  two  possible  solutions,  specular 
in  each  element.  This  is  not  in  contradiction  with  the  uniqueness  result  for  Generalized  Gauss 
quadrature  rules  as  available  in  [12],  since  the  functions  that  we  are  integrating  do  not  satisfy 
the  hypothesis  of  this  result,  in  particular  they  are  not  a  Chebychev  System.  On  the  other  hand, 
for  the  case  (m  —  q)  odd,  we  have  several  admissible  solutions.  Among  these,  there  is  one  that 
is  symmetric  in  the  element  and  our  choice  is  to  select  exactly  this  one.  To  this  purpose,  we 
complement  (3.5)  with  the  condition  x[  +  x! ,  =  1/2. 

We  introduce  now  the  functions  {pj}j= i,...,m-q  such  that 


r]j(x)  =  fij'\x  +  i) 


Tjj  (x)  =  y,l+1\x  +  i)  +  ip%’l+1,(x  +  i  +  1) 


(M+i ), 


j  =  1, . . . ,  m  -  2q  -  1 


j  =  m  —  2g, . . . ,  m  —  q 
and  h  —  j  —  m  +  2g  +  1, 


(3.6) 


for  all  x  G  [0, 1].  The  functions  r/j  are  (m  —  q )  functions,  independent  of  i,  with  support  on  [0, 1] 
such  that: 

f\  t 


'o 


rjj(x)dx  =  j  y  (x)dx  j  =  1, . . . ,  m  —  2g  —  1 


-,«( 


Ji-1 

r  1  pi- 1-1 

/  r]j(x)dx  =  /  (x)dx  j  =  m  —  2g, . . . ,  m  —  q ;  h  =  j  —  m  +  2g  +  1 

'0  Ji- 1 


10 


Figure  3:  Functions  to  be  exactly  integrated  by  the  quadrature  rule  in  Problem  3.2,  namely  functions  rjj(x)  of 
Equations  (3.6).  A  choice  of  the  resulting  quadrature  points  is  plotted  on  the  x  axis.  On  the  left  is  the  case  with 
in  =  10,  k  =  3,  q  =  4,  on  the  right  is  the  case  with  m  =  10,  k  =  3,  q  =  3. 


Figure  3  shows  for  two  choices  of  m,k  and  q  the  functions  rjj(x)  and  the  resulting  quadrature 
points.  Setting 

Mj  =  /  r)j(x)dx , 

Jo 

we  can  rewrite  Problem  3.1  in  an  equivalent  form  as  follows 
Problem  3.2.  Find  xj  E  [0, 1]  and  wf  >  0  ,  l  —  1, . . . ,  n1  such  that 

n1 

^2wiVjixi)  =Mj. 

i=i 

□ 


This  is  the  only  non-linear  problem  that  needs  to  be  solved  in  our  final  algorithm.  We  stress 
that  the  dimension  of  this  non-linear  problem  only  depends  on  the  degree  m  and  regularity  q 
and  it  is  independent  of  k,  i.e. ,  it  is  independent  of  the  number  of  elements  introduced  in  the 
discretization. 

However,  when  m  is  large,  Problem  3.2  can  be  ill-conditioned,  and  for  this  reason  we  deal 
with  it  following  the  approach  proposed  in  [12].  Assuming  for  simplicity  (to  —  q)  to  be  even1,  we 

1  Remember  that  in  the  case  (to  —  q)  odd  we  consider  one  new  equation  imposing  symmetry  and  recover  the 
square  dimension.  This  modifies  the  functional  <G  described  below,  although  the  Jacobian  can  be  calculated  also  in 
this  case. 
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introduce  the  following  functional  €  defined  as 


€  :  M2"7  ->  R2nI  ,  2 n1  =  m-q 


where 

n1 

CiyCr,.  ...,xni,w  i, . . . ,  wni)  =  ^2  wiVj(xi)  ~  Mj. 

i=i 

With  this  notation,  the  solution  of  Problem  3.2  is  a  zero  of  the  above  functional: 


€(x7,  w7)  =  0. 


Now,  the  Jacobian  of  functional  €  can  be  easily  constructed  explicitly  since  the  derivatives  of  the 
basis  functions  are  known,  in  fact: 

dGj  ,  OQj 

Thus,  Newton  iterations  can  be  constructed,  ft  is  important  to  observe  that  the  choice  of  the 
starting  point  is  crucial  to  obtain  convergence,  due  to  the  possible  ill-conditioning  of  the  Jacobian 
matrix.  As  proposed  in  [12],  we  use  a  continuation  algorithm  to  obtain  better  conditioned  sub¬ 
problems.  Therefore,  we  are  led  to  solve  the  following: 


(on 


and  a  sequence  of  increasing  an  €  (0,1],  it  = 


Problem  3.3.  Fix  , . . . ,  x 1  /  ,w{  , 

1, . .  with  anit  =  1. 

For  all  it  =  1, ... ,  nit,  find  (xf \  . . . ,  x^\w^\  . . . ,  w$)  s.t.  <G(tt)(4rt), . . . ,  . . . ,  u/f ) 

0,  where  we  define: 


■  ■  -,xni,w  i, . .  .,wni) 


n1 

y^/wirijjxi)  -  ait 

/  T]j{x)dx 

-  (1  -  <%) 

n1 

^2lwltt~1)r}j(x\%t~1)) 

i=i 

J  0 

i=i 

We  propose  to  use  as  initial  choice  (x^\  . . . ,  . . . ,  w®})  the  Gauss  nodes  and  equal 

weights,  and  a  vector  of  parameters  ait  fitted  near  zero. 

Summarizing,  the  solution  of  Problem  3.1  is  given  by  xj  =  x[nxt\  and  wj ’*  =  w[nu\  for  l  = 
l,...,n7  =  (m-q)/ 2. 
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Figure  4:  Functions  to  be  integrated  by  the  boundary  quadrature  rule,  namely  functions  (x)  and 

'  (x)  truncated  to  the  first  element.  Note  that  the  rule  has  to  give  the  inexact  result  for  the  latter,  see  Equation 
(3.7).  On  the  left  is  the  case  with  to  =  10,  k  =  3,  q  =  4,  on  the  right  is  the  case  with  to  =  10,  k  =  3,  q  =  3. 


3.3.  Quadrature  on  an  open  uniform  knot  vector 

As  commented  in  Section  3.1,  the  spline  space  S™  generated  by  the  open  uniform  knot  vector 
is  larger  than  the  space  S™  considered  in  the  previous  section,  since  it  contains  also  the  boundary 
functions  and  . .  . ,  V’q+V  Then,  the  quadrature  rule  of  Section  3.2  needs  to  be 

modified  properly  at  the  boundary. 

For  exposition  purposes,  we  focus  on  the  left  boundary;  thus,  we  have  to  construct  a  quadrature 
rule  for  the  basis  functions  whose  support  contains  [0,1].  Following  the  notation  introduced  in 
Section  3.1,  such  basis  functions  coincide  with  the  first  m  +  1  basis  functions  of  the  spline  space 
S™,  that  is 


{A}i=l,.,m+l  =  {^l\  •  •  •  ,  tpq+l,  •  •  •  ,  ■  ■  ■  ,  ^q+l}  ■ 


d!)  „/,(!) 


(1) 


,(1,2) 


',(1,2)  1 


and  these  functions  are  plotted  in  Figure  6  for  two  choices  of  m,  k  and  q. 

We  note  that  that  the  set  {'0i|[o,i]}i=i,...,m+i  is  a  basis  for  the  polynomials  of  degree  m  on  the 
first  element  [0,1]. 

Our  aim  is  to  solve  the  following: 

Problem  3.4  (Boundary  Quadrature).  Find  nB  quadrature  points  xf  G  [0,1],  and  weights  wf, 
l  =  1, . . . ,  nB  such  that: 

nB  nI 

y ^wfipj(xf)=  I  ipj(x)dx y  wl’^i^xl'2),  i  =  l,...,m  +  l,  (3.7) 

i=i  i=i 
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Figure  5:  The  resulting  quadrature  on  the  boundary:  the  function  g(x)  defined  in  (3.8)  and  the  quadrature  rule 
in  the  case  of  m  +  1  =  11  Gauss  nodes.  On  the  left  in  the  case  m  =  10;  k  =  3;  q  =  4  and  on  the  right  in  the  case 
m  =  10;  k  =  3;  q  =  3. 


where  xl’2,wj’2  are  the  quadrature  points  and  weights  from  Problem  (3.1). 


□ 


In  order  to  represent  the  right  hand  side  of  (3.7)  in  a  more  convenient  way,  we  introduce  a 
function  g{x)  such  that 

r  1  r2  n1 


ipi(x)g(x)dx  =  /  ipj(x)dx  -  W  irj  V 


(3.8) 


1=1 


JO  Jo 

Function  g{x)  can  be  chosen  as  an  m-degree  polynomial  and  can  be  easily  computed,  see  Section 
3.4.  Then  we  have  the  following  weighted  quadrature  problem  equivalent  to  Problem  3.4. 

Problem  3.5.  Find  nB  quadrature  points  x f  G  [0, 1],  and  weights  wf ,  l  =  1, . . .  ,nB  such  that: 

nB 


=  /  ifj(x)g(x)dx  Vj  =  1, . . .  ,m  +  1. 


i=i 


(3.9) 

□ 


Solving  Problem  3.5  is  equivalent  to  finding  a  quadrature  rule  with  respect  to  the  weighted 
measure  g(x)dx,  requiring  exactness  for  polynomials  of  degree  m  on  [0,1].  Since  g(x)  has  in 
general  oscillating  sign,  the  existence  of  a  quadrature  rule  with  [m  +  l)/2  points  is  not  guaranteed 
in  this  context.  However,  we  can  use  an  (m  +  l)-nodes  standard  Gauss  rule  {xf  ,wf ):  it  integrates 
exactly  the  polynomials  of  degree  2 m  +  1,  and  in  particular,  all  the  products  fjJ(x)g(x).  Thus  we 
have: 

1  m+l 

/  ifj(x)g(x)dx  =  ^wf^j{xf)g{xf)  .  (3.10) 


i=i 
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Problem  3.5  is  then  solved  for  nB  =  m  +  1  by  setting  wf  =  wfg(xf )  and  xf  =  xf . 

In  Figure  5  the  quadrature  rule  is  plotted.  Notice  that  some  of  the  final  weights  are  negative. 
Related  stability  issues  are  discussed  in  Remark  3.1. 

The  advantage  of  this  procedure  is  that  g(x)  can  be  accurately  and  easily  computed.  For  this 
purpose,  we  introduce  an  auxiliary  basis  {0}  for  the  space  of  polynomials  of  degree  m  on  [0, 1]. 
Choosing  Q  =  ipj  our  algorithm  (see  Section  3.4)  works  up  to  machine  precision  for  m  <  16.  If  we 
choose  instead  Cj(x)  as  the  Legendre  polynomials  we  can  go  up  to  degree  m  —  32. 

To  summarize,  suppose  to  have  solved  Problem  3.2  and  Problem  3.4  (for  the  left  boundary 
(xfB,wfB)  =  (xf,wf)  and  for  the  right  boundary  ( xfB,wfB )  in  the  analogous  way),  then  the 
quadrature  rule  for  exact  integration  of  all  /  e  S™  is 

nB  k— 2  /  n1  \  nB 

nx)  <BftfB)  +  E  E  >»//(*? + o  +  E  <Bn*RB)- 

i=i  i= i  y  i=i  J  i=i 

We  remark  that  if  m  —  q  is  odd,  the  two  boundary  problems  are  specular  and  (xBB ,  wBB)  are 
immediately  obtained  from  (xBB  ,wBB). 

Remark  3.1.  Stability  of  a  quadrature  rule,  and  its  convergence  on  C°  integrand  functions  as  the 
number  of  quadrature  points  increases,  is  guaranteed  when  the  sum  of  the  weight  absolute  values  is 
uniformly  bounded,  see  for  example  [11].  In  the  case  of  (3.3),  we  easily  have  that 

n1 

Yd\u,l\  = 

1=1 

since  internal  weights  w]  are  positive  ( recall  Problem  3.1).  Furthermore,  in  all  our  numerical  tests, 
up  to  m  =  32,  we  have  always  observed  that  Ym=i  \wLB\  is  bounded,  therefore,  despite  that  some 
quadrature  weights  on  the  boundary  elements  may  be  negative,  the  stability  of  the  overall  procedure 
is  guaranteed. 

Remark  3.2.  One  can  also  consider  alternative  procedures  to  avoid  negative  weights.  As  in  Prob¬ 
lem  3.2,  one  can  try  to  numerically  compute  the  quadrature  points  xf  and  weights  wf  in  Problem 
3.5,  constraining  the  wf  to  be  positive.  Existence  of  the  solution  is  not  guaranteed  but  our  experi¬ 
ence  shows  that  a  solution  can  be  found  with,  e.g.,  nB  —  m  +  1,  as  in  (3.10). 

3-4-  Quadrature  algorithm 

A  pseudo-code  for  the  numerical  evaluation  of  the  quadrature  points  and  weights  is  proposed 
in  this  section,  again  for  simplicity  in  the  case  of  m  —  q  even.  The  output  is  the  nodes  and  weights 
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for  internal,  left  boundary  and  right  boundary.  All  the  three  sets  of  nodes  are  calculated  in  [0, 1], 
thus,  once  computed,  they  have  to  be  properly  shifted  in  the  patch. 


Input:  m  =  spline  degree;  q  =  spline  regularity; 

Set  n1  =  (m  —  q)/ 2,  nB  —  m  +  1; 

{Computation  of  the  quadrature  rule  on  internal  elements} 

Set  a  be  a  vector  of  nit  ordered  real  numbers  s.t.  a(i)  >  0 ,  a(nit )  =  1. 

Calculate  A4j,  exact  integrals  of  functions  i]j(x)  defined  in  (3.6); 

Initialize  xp ,  w j0'1; 

for  i  —  1, . . . ,  nit  do  {see  Problem  3.3} 

Compute  wp ,  xp  such  that 

wP  >  0,  xp  G  [0, 1],  V/  =  1, . . . ,  n1  and 
JP  wPrjj(xP)  -  a(i)Mj  -  (1  -  a(i))  w^p^x^)  =  0. 

i=i  i=i 

end  for 

Set  (xf,wf)  =  (xpit]  .wp^). 

{Computation  of  the  quadrature  rule  on  the  left  boundary  element}. 

Calculate  A  =  (aij)ij=l!_tm+1 ,  a{j  =  f*  ilii(x)(j(x)dx; 

Calculate  A*  =  JQ2  ipi(x)dx  —  X]T=i  ^^(xf  +  1); 

Solve  the  linear  problem  A/3  =  A  and  set  g{x )  =  Ypj=i  PjCj(x)> 

Calculate  a  n^-points  quadrature  rule  =  (xf,wf)  of  degree  of  exactness 

m  for  the  weighted  measure  g(x)dx.  {e.g.,  use  (3.10)}. 

{Computation  of  the  quadrature  rule  on  the  right  boundary  element}. 

Calculate  in  a  similar  way  (xfB  ,wPB). 

return  Internal  n7-points  quadrature  rule  (xf,tcf)  and  two  boundary  nB-points 
quadrature  rules  (xBB ,wBB),  (xfB ,wBB). 

In  all  our  tests  system  (3.11)  is  solved  by  the  MATLAB  Optimization  Toolbox  fsolve  routine 
(see  [26]),  setting  a  machine-precision  tolerance  on  the  residual  norm.  The  complete  MATLAB 
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implementation  for  the  calculation  of  the  interior  and  boundary  quadrature  is  available  in  the 
Contributions  section  at  geopdes.sourceforge.net. 

Below,  we  also  report  the  pseudo  code  for  the  alternative  calculation  of  the  quadrature  rule  at 
the  boundary,  with  enforcement  of  positivity  of  the  weights,  as  described  in  Remark  3.2. 


Set  nB ,  e.g.,  nB  —  m  +  1; 

{Alternative  computation  of  the  quadrature  rule  on  the  left  boundary  element} 
Calculate  A*  =  JQ2 f>i(x)dx  —  Y^Ti=iwi/lPi(xi  +  1)> 

Compute  wi,  Xi  such  that 
' 

u>i  >  0,  xi  e  [0, 1],  V/  =  1, . . . ,  nB  and 

<  nB  (3.12) 

E  -  Ai  =  0,  Vi  =  1, . . . ,  m  +  1 

,  i= i 

Set  (xfB,wBB)  =  (. xhwt ). 


3.5.  Optimality  of  the  new  quadrature  rule 

In  the  previous  sections  we  have  proposed  a  new  quadrature  rule,  which  can  be  computed  in 
an  efficient  and  stable  manner.  We  now  want  to  discuss  the  optimality  of  the  proposed  quadrature 
rule. 

We  observe  the  following. 


•  In  each  space  dimension,  the  proposed  rule  uses  ( k  —  2)  — - — 
points,  with  m  =  2 p  and  q  —  r  —  2,  when  NURBS  of  degree  p  anc 
for  the  IGA  of  an  elliptic  second  order  problem. 


+  2(m  +  1)  quadrature 
regularity  r  are  selected 


•  The  optimal  quadrature  rule  proposed  in  [21]  would  require 
points  in  each  dimension. 


k 


m  —  q 
2 


quadrature 


•  Element-wise  Gauss  quadrature  needs  k 


m  +  l 
2 


points  per  dimension  instead. 


According  to  these  considerations,  in  Figure  6  we  plot  the  number  of  quadrature  points  (i.e.,  the 
computational  cost)  for  the  three  different  quadrature  rules  reported  above  in  the  case  of  a  uniform 
three-dimensional  mesh  of  k  x  k  x  k  elements.  The  figure  shows  that  the  quadrature  rule  proposed 
in  this  paper  has  approximately  the  same  number  of  quadrature  points  as  the  optimal  one  from 
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P  k 


Figure  6:  Comparison  of  the  number  of  quadrature  points  needed  for  various  quadrature  rules  on  a  (k  x  k  x  k)- 
element  three-dimensional  mesh.  On  the  left:  p  =  1, . . . ,  16  (to  =  2 p),  q  =  p  —  2,  and  k  =  4, . . . ,  100,  on  the  right: 
p  =  6  (to  =  12),  q  =  p  —  2  =  4,  and  k  =  4, . . . ,  100. 

[21],  which  is  approximately  1/8  the  number  of  quadrature  points  of  the  element-wise  Gauss  rule, 
for  large  k  and  p. 

However,  we  recall  that  the  optimal  quadrature  rule  from  [21]  is  in  practice  impossible  to 
obtained  for  the  cases  considered  in  Figure  3.5  (i.e.  p  up  to  16  and  k  up  to  100),  while  the  com¬ 
putation  of  quadrature  points  and  weights  for  the  new  rule  is  straightforward  with  the  algorithm 
discussed  in  Section  3.4. 

4.  Applications 

We  now  present  some  applications  of  the  proposed  new  quadrature  rule.  In  particular,  we  prove 
its  efficiency  versus  the  element-wise  Gauss  quadrature  rule,  for  boundary  value  problems  in  two- 
dimensions  solved  with  a  standard  IGA  formulation  (see  [13]  for  details).  We  remark  that  2D  (and 
3D)  integration  rules  are  simply  constructed  as  tensor  products  of  ID  rules.  We  moreover  discuss 
how  to  solve  examples  where  some  types  of  non-uniform  meshes  are  considered.  In  particular,  we 
show  the  behavior  of  the  new  rules  on  an  example  where  a  non-uniform,  anisotropically  graded 
mesh  is  used  to  capture  a  solution  exhibiting  thin  layers. 

4-1.  Numerical  solution  of  a  Poisson  problem  on  a  quarter  of  an  annulus 

As  a  first  example,  a  problem  involving  a  non-trivial  geometry  that  can  be  exactly  represented 
by  NURBS  is  considered.  With  this  aim,  we  focus  on  a  Poisson  problem  defined  on  a  quarter  of 
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an  annulus  fl  (see  Figure  7)  with  an  internal  radius  R\  =  1  and  an  external  radius  R2 
we  solve  the  model  problem: 

j  —A u  =  /,  \/x  G  hi, 

[  u\on  =  0, 

with 


4,  and 

(4.1) 


/  =  (2a:4  —  50a:2  —  50 y2  +  2y4  +  4a:2a/2  +  100)  sin(x)  sin (y)  + 

+  (68a:  —  8a:3  —  8  xy2)  cos(a:)  sin(y)  +  (68  y  —  8  y3  —  8  yx2)  cos  (y)  sin  (a;), 

such  that  the  exact  solution  is 

u  =  (x2  +  y2  —  l)(a:2  +  y2  —  16)  sin(x)  sin(y). 


Figure  7:  Poisson  problem  on  a  quarter  of  an  annulus:  geometry  of  the  domain  Q. 

We  find  an  approximation  of  the  solution  by  a  standard  IGA  Galerkin  formulation,  using 
for  integration  the  proposed  new  quadrature  rules,  as  well  as,  for  comparison  reasons,  a  standard 
element-wise  Gauss  quadrature.  The  problem  is  solved  for  basis  function  degrees  p  =  q  =  2, . . . ,  5  in 
both  parametric  directions  (and  maximum  inter-element  regularity),  using  control  nets  consisting 
of  25  x  25  up  to  100  x  100  control  points.  In  Figure  8,  we  show  the  convergence  curves  for  the 
L2-norm  of  the  relative  error  with  respect  to  the  exact  solution.  The  theoretically  expected  rates 
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of  convergence  (i.e. ,  p  +  1)  are  obtained.  In  the  same  Figure,  we  also  report  the  error  computed 
with  respect  to  the  numerical  solution  obtained  using  element-wise  Gauss  quadrature.  It  can  be 
seen  that  this  error  is  always  negligible  with  respect  to  the  approximation  error  (a  saturation  is 
observed  for  high  orders,  when  this  error  is  in  the  range  of  machine-precision). 


— e—  p  =  2  (vs  exact) 
— e—  p  =  3  (vs  exact) 
— e— p  =  4  (vs  exact) 
— e— p  =  5  (vs  exact) 
— e-p  =  2  (vs  Gauss) 
-^^p  =  3  (vs  Gauss) 
MH-p  =  4  (vs  Gauss) 
-^^p  =  5  (vs  Gauss) 


Figure  8:  Poisson  problem  on  a  quarter  of  an  annulus:  convergence  curves  for  the  L2-norm  of  the  relative  error 
(j)  =  q  =  2, . . . ,  5),  using  for  integration  the  proposed  new  quadrature  rules;  errors  are  computed  both  with  respect 
to  the  exact  solution  (thick  curves)  and  to  the  numerical  solution  obtained  using  element-wise  Gauss  quadrature 
(thin  curves). 


In  Figure  9  we  report  the  total  number  of  points  used  for  integration  by  the  two  quadrature 
strategies.  In  particular,  on  the  top,  we  show  the  number  of  quadrature  points  as  a  function  of  the 
number  of  control  points  for  a  fixed  basis  function  degree  p  —  q  —  4  (and  maximum  inter-element 
regularity)  in  both  parametric  directions,  while,  on  the  bottom,  we  show  the  number  of  quadrature 
points  as  a  function  of  the  basis  function  degrees  p  =  q  (always  considering  maximum  regularity) 
for  a  fixed  100  x  100  control  net.  The  advantage  of  the  new  rules  is  clear,  in  particular  for  high 
degrees  (implying  high  inter-clement  regularity)  and  when  the  number  of  control  points  is  not  too 
low,  that  is,  when  the  number  of  boundary  elements  (i.e.,  the  number  of  elements  where  the  rules 
are  not  optimal)  is  not  large  with  respect  to  the  number  of  internal  elements. 
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Figure  9:  Poisson  problem  on  a  quarter  of  an  annulus:  number  of  quadrature  points  as  a  function  of  the  number  of 
control  points  for  a  fixed  basis  function  degree  p  =  q  =  4  in  both  parametric  directions  (top) ;  number  of  quadrature 
points  as  a  function  of  the  basis  function  degrees  p  =  q  for  a  fixed  100  x  100  control  net  (bottom). 
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Finally,  in  Figure  10,  we  show  the  distribution  of  the  points  required  by  the  two  quadrature 
strategies  in  the  case  p  =  q  =  4  when  a  25  x  25  control  net  is  adopted.  The  advantage  of  using  the 
new  rules  is  clear  (5625  points  versus  the  11025  needed  by  Gauss  quadrature),  despite  this  case 
not  being  a  very  fine  mesh  (where,  as  seen  above,  such  an  advantage  would  be  amplified). 

4-2.  Extension  to  non-uniform  meshes 

The  theoretical  part  of  the  present  paper  focuses  only  on  the  construction  of  efficient  quadrature 
rules  for  the  case  of  uniform  knot  vectors.  In  this  section,  we  sketch  how  to  extend  the  proposed 
new  quadrature  rules  to  some  particular,  but  useful,  cases  of  non-uniform  knot  vectors.  We  however 
remark  that,  despite  the  fact  that  many  of  the  typical  cases  of  mesh  refinement  in  IGA  can  be 
covered,  as  discussed  in  the  following,  a  solution  for  more  general  non-uniform  situations  is  not 
within  the  scope  of  this  work. 

4-2.1.  Geometrically  graded  mesh 

The  problem  of  the  extension  of  the  proposed  quadrature  formula  to  non-uniform  knot  vectors 
is  made  difficult  by  the  issue  of  the  proper  integration  of  transmission  functions,  which  relies 
strongly  on  their  translation-invariant  property.  The  removal  of  the  hypothesis  of  knot  uniformity 
implies  that,  in  general,  such  a  translation- invariant  property  no  longer  holds.  However,  since  all 
the  functions  to  be  exactly  integrated  on  interior  elements  are  supported  over  only  two  elements,  a 
way  to  simply  construct  the  new  quadrature  rules  for  a  non-uniform  knot  vector  is  to  constrain  the 
ratio  between  the  lengths  of  two  consecutive  knot  spans  to  be  fixed  (i.e.,  a  “geometrically  graded” 
knot  vector  is  considered).  In  this  case,  when  considering  pairs  of  consecutive  elements,  it  is  clear 
that  transmission  functions  always  have  the  same  structure  and,  as  a  consequence,  the  previously 
described  procedure  for  computing  quadrature  points  and  weights  can  be  simply  extended  to  this 
situation. 

This  is  indeed  an  interesting  case  from  the  point  of  view  of  applications,  since  geometrically 
graded  meshes  are  among  the  most  used  categories  of  non-uniform  meshes  and  are,  in  particular, 
important  when  solution  layers  have  to  be  captured  (see,  e.g.,  the  example  in  the  next  section). 

4-2.2.  Quadrature  multi-patch 

We  define  quadrature  multi-patch  as  the  case  of  a  mesh  made  of  “sub-patches”  where  a  dif¬ 
ferent  (but  unique  within  each  single  sub-patch)  quadrature  rule  is  used.  The  sub-patches  are 
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Figure  10:  Poisson  problem  on  a  quarter  of  an  annulus:  distribution  of  integration  points,  in  the  case  of  basis 
function  degree  p  =  q  =  4  in  both  parametric  directions  and  a  25  x  25  control  net,  for  the  proposed  new  quadrature 
(top  left)  and  the  element-wise  Gauss  quadrature  (top  right).  On  the  bottom,  details  of  the  bottom  right  corner  of 
the  plots  above  are  reported. 
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interconnected  by  the  so-called  “transition”  zones.  Quadrature  multi-patch  can  be  easily  treated, 
as  well. 

In  this  situation,  the  suggested  solution  is  to  perform  integration  considering  each  sub-patch 
individually,  as  if  it  was  not  connected  to  the  others.  Thus  the  proper  internal  quadrature  rule 
has  to  be  used  for  the  internal  elements  of  each  sub-patch,  while  boundary  integration  has  to  be 
performed  for  the  elements  lying  on  the  mesh  borders  as  well  as  for  transition  zones.  An  example 
is  given  in  Figure  11,  where  a  mesh  constituted  by  6  different  sub-patches  is  constructed  as  the 
tensor  product  of  piecewise  uniform  and  geometrically  graded  knot  spans.  The  elements  where 
boundary  quadrature  is  in  order  are  those  adjacent  to  the  black  thick  lines. 

In  particular,  this  possibility  allows  to  perform  integration  on  meshes  that  are  composed  of 
uniform  zones  (even  of  different  mesh-sizes)  and  of  non-uniform,  geometrically  graded  zones  (as 
described  in  the  previous  section),  i.e. ,  all  the  possible  combinations  of  the  quadrature  rules  consid¬ 
ered  in  the  present  paper  are  possible.  As  a  consequence,  integration  can  be  efficiently  performed 
on  a  large  variety  of  mesh  situations,  suitable  to  solve  many  of  the  problems  typically  tackled 
by  means  of  tensor  product  NURBS-based  IGA.  An  interesting  example  is  reported  in  the  next 
section. 

We  remark  that  this  procedure  is  cost-effective  only  if  the  number  of  transition  elements  is 
significantly  less  than  internal  elements.  If  this  is  not  the  case,  the  macro-element  quadrature 
presented  in  [21]  is  preferred. 

4-3.  Numerical  solution  of  a  reaction- diffusion  problem  with  layers,  on  a  suitable  non-uniform 
mesh 

In  this  section  we  make  use  of  the  non-uniform  quadrature  strategies  discussed  in  the  previous 
section,  to  compute  the  integrals  necessary  to  solve  the  following  reaction-diffusion  problem  defined 
on  the  bi-unit  square  Id  =  [0,  l]2: 

{du 

_io-vAt,  +  — =i,  v*  e  a  =  [o,  i]2,  (42) 

u\on  =  0. 

The  exact  solution  of  such  a  problem  is  a  ramp  of  unit  slope  along  the  x-axis,  showing  two  layers 
at  y  =  0  and  y  —  1,  and  a  third,  sharper  layer  at  x  —  1. 

In  order  to  solve  this  problem  with  a  uniform  mesh,  a  very  small  mesh  size  has  to  be  selected 
in  order  to  correctly  capture  the  layers.  Therefore,  a  non-uniform  mesh,  graded  towards  the 
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Figure  11:  Quadrature  multi-patch:  example  of  a  mesh  consisting  of  6  different  sub-patches,  constructed  as  the 
tensor  product  of  piecewise  uniform  and  geometrically  graded  knot  spans.  The  elements  where  boundary  quadrature 
is  necessary  are  those  adjacent  to  the  black  thick  lines. 

layers,  is  certainly  preferred.  Using  the  non-uniform  quadrature  strategies  previously  introduced, 
integration  over  a  mesh  like  the  one  reported  in  Figure  12  is  possible,  allowing  capturing  the  layers 
of  the  solution  with  a  reasonable  number  of  elements.  The  solution  computed  with  a  standard 
Galerkin  IGA  formulation  (degree:  p  =  q  =  3,  control  net:  21  x  21)  on  such  a  mesh,  using  the  new 
quadrature  rules,  is  reported  in  Figure  12  as  a  color  map,  while  in  Figure  13  a  3D  plot  illustrates 
more  clearly  how  the  layers  are  reproduced  by  the  numerical  solution. 

5.  Summary  and  Conclusions 

We  have  considered  the  development  of  efficient  quadrature  rules  for  NURBS-based  Isogeomet¬ 
ric  Analysis.  We  have  focused  on  rules  for  arrays  that  frequently  arise  in  Finite  Element  Analysis, 
namely,  mass,  stiffness  and  advection  matrices.  We  have  reduced  the  quadrature  problem  to  the 
exact  integration  of  certain  spaces  of  B-Splines  in  one  dimension  that  have  essential  translation 
invariance  and  localization  properties,  and  for  this  cases  we  have  developed  accurate  rules  for  pe¬ 
riodic  B-Splincs  and  B-Splines  constructed  for  open  knot  vectors.  Both  uniform  and  geometrically 
scaled  knot  spacing  have  been  considered.  We  have  presented  a  simple  and  efficient  algorithm  for 
constructing  the  rules  and  we  have  illustrated  its  use  on  boundary  value  problems,  including  one 
whose  solution  possesses  sharp  boundary  layers.  The  way  to  use  the  rule  on  domain  consisting  of 
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Figure  12:  Reaction-diffusion  problem  with  layers:  adopted  non-uniform  mesh  and  color  map  of  the  corresponding 
numerical  solution  (degree:  p  =  q  =  3,  control  net:  21  x  21). 
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Figure  13:  Reaction-diffusion  problem  with  layers:  3D  plot  of  the  computed  numerical  solution  (degrees:  p  =  q  =  3, 
control  net:  21  x  21). 

multiple  patches  has  also  been  described.  Our  results  indicate  that  the  rules  compare  very  favor¬ 
ably  with  optimal  rules,  which  are  at  the  very  best  difficult  to  obtain  and  in  some  cases  almost 
impossible.  The  new  rules  are  of  course  more  efficient  than  Gauss  quadrature  repeated  on  knot 
spans.  We  believe  this  work  takes  another  step  forward  in  the  development  of  specialized  quadra¬ 
ture  rules  for  NURBS-based  Isogeometric  Analysis.  An  obvious  challenge  is  to  develop  rules  for 
general  non-uniform  knot  spacing.  We  plan  to  pursue  this  in  future  work. 
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